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<N Abstract 

Q We report a study of the homogeneous isotropic Boltzmann equation for an open system. We 

u seek for nonequihbrium steady solutions in presence of forcing and dissipation. Using the language 

d 

•jH of weak turbulence theory, we analyze the possibility to observe Kolmogorov-Zakharov steady dis- 

I— I tributions. We derive a differential approximation model and we find that the expected nonequilib- 

rium steady solutions have always the form of warm cascades. We propose an analytical prediction 
for relation between the forcing and dissipation and the thermodynamic quantities of the system. 
Specifically, we find that the temperature of the system is independent of the forcing amplitude 
and determined only by the forcing and dissipation scales. Finally, we perform direct numerical 
simulations of the Boltzmann equation finding consistent results with our theoretical predictions. 
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I. INTRODUCTION 



Systems in a steady state are characterized by observables that do not change in time; 
they can be either in equihbrium or out of equihbrium. Systems in nonequihbrium steady 
states have net currents (fluxes): examples of nonequihbrium steady-state systems include 
an object in contact with two thermal sources at different temperatures, for which the cur- 
rent is a heat flux; a resistor with electric current flowing across it; the kinesin- microtubule 
system, for which kinesin motion is the current. Most biological systems, including molecu- 
lar machines and even whole cells, are in nonequihbrium states [IJ. In particular, biological 
systems rely on a continuous flux of energy and/or particles supplied by some proper envi- 
ronmental reservoirs. 

In statistical mechanics, investigating the general properties of a system in contact with 
reservoirs, namely an open system, is a long lasting problem (e.g. see the second problem 
discussed by E.H. Lieb on the occasion of the award of the Boltzmann medal [2]), even 
though these theoretical challenges are sometimes neglected in applied engineering at large. 
The difficulties arise from the fact that finding the large deviation functional for a stationary 
state with fluxes is still an open problem (see [3] and references therein). In the present work, 
for focusing our attention and considering an affordable goal, we consider the kinetic theory 
of gases. In particular, we consider a system composed of a large number of interacting 
particles, comparable to the Avogadro number. The Boltzmann kinetic equation (BKE) 
describes the time evolution of the single-particle distribution function, which provides a 
statistical description of the positions and velocities (momenta) of the gas molecules. This 
integro-differential kinetic equation, proposed by Boltzmann at the end of the XIX century, 
has been derived starting from the phase-space Liouville equation, assuming the stosszahl 
ansatz [1]. Its equilibrium state, which maximizes the entropy measure, is the Maxwell- 
Boltzmann distribution. In case of small deviations from the local equilibrium, it is possible 
to systematically derive hydrodynamic equations for macroscopic quantities of the system; 
e.g., in the lowest order approximation for small departures from equilibrium, the Navier- 
Stokes equations 

Kinetic equations have also been studied in the framework of wave turbulence theory 
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[5] where it has been shown that other solutions with respect to thermodynamic solutions 
can be stationary states of the system, in case of external forcing and dissipation. These 
distributions, which have usually the form of power-laws in momentum space, are called 
Kolmogorov-Zakharov (KZ) and they represent constant flux of conserved quantities similar 
to the Kolmogorov energy cascade in strong Navier- Stokes turbulence P, [7j. These solu- 
tions, named cascade solutions, become important when considering an open system, i.e. 
with forcing and dissipation terms. They have been studied for a great variety of weakly 
nonlinear dispersive models: examples can be found in water waves [SHTU]. internal waves 
[TT] . nonlinear optics [12], Bose- Einstein condensation [T3HI5], magnetohydrodynamics [TB] . 

An out of equilibrium description of the Boltzmann equation using the KZ solutions was 
first devised in [17] considering different types of interaction potential between particles. 
Problems of interaction locality scale-by-scale and wrong flux direction were pointed out. In 
particular in [18] Kats showed that for all realistic physical situations the direction of the 
cascades in the system is always in the wrong orientation with respect to the one predicted 
by the Fj0rtoft theorem [23] • When a formal KZ solution has a flux direction contradicting 
with the Fj0rtoft theorem, this spectrum (even if local) cannot be established because it 
cannot be matched to any physical forcing and dissipation at the ends of the inertial range. 
For example in [12] , the particle cascade KZ solution was found to be of this type in the two- 
dimensional nonlinear Schrodinger equation model the authors argued that in this case the 
KZ solution is not achievable and a mixed state, with both a cascade and a thermodynamic 
components were proposed. Another example of mixed cascade-thermodynamic states can 
be found in the context of three-dimensional Navier-Stokes turbulence [19] , where such mixed 
states were called warm cascades [M] . 

The present manuscript will focus on warm cascades found in the homogenous isotropic 
Boltzmann equation (HIBE) and in particular it will answer to the following important 
questions. 

• What is precisely the relation between the conserved quantity fluxes and the thermo- 
dynamics quantities of the system? 

• How does this relation depends on the forcing and dissipation rates and acting scales? 

To answer the above questions we will perform numerical simulations of the homogeneous 
isotropic Boltzmann equation with forcing and dissipation. We will then use a diffusion 
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approximation model (DAM) to derive analytical predictions on how the thermodynamic 
quantities, temperature and chemical potential, are related to fluxes, forcing and dissipative 
scales. We will then test these predictions by numerically simulating both DAM and the 
complete homogenous isotropic Boltzmann equation. 

The work is organized as follows: in Section ITT] we review the properties of the Boltzmann 



equation for the homogeneous isotropic case; in Section III we introduce DAM and we derive 



the analytical predictions; Section [IV] is dedicated to numerical results of DAM and HIRE; 
in Section |V] we draw the conclusions. A set of Appendixes also provide detailed calculations 
of those results which are briefly reported in the main text. 



II. THE BOLTZMANN KINETIC EQUATION 

The Boltzmann kinetic equation describes the time evolution of the single-particle dis- 
tribution function, which provides a statistical description for the positions and momenta 
of the gas molecules: the function n(x, k, t) express a probability density function in the 
one-particle phase space x with respect to time, where d is the dimension. Note that 
we denote the momentum variable with the letter k instead of the conventional p to follow 
the common notation of wave turbulence [5]. The Boltzmann equation takes the following 
form: 

— (x,ki,)f:) H — (x,ki,t) = /co//(x,ki,t), (1) 

ot m ox 

where 

/ + 00 
[^(X) kg, t)n(x, k4, t) - n(x, ki, t)n{x, ka, t)] (ik2c/k3c/k4 (2) 
-oo 

sums the effect of the two-body collisions of particles with all possible values of momenta. 
The form of the collision integral we are reporting is equivalent to the standard one and 
corresponds to Eq. (4.18), page 64 in Cercignani's book Here W describes synthetically 
the scattering amplitude transition 2 — )■ 2 as a function of the momenta of the interacting 
particles. As we consider elastic collisions, the general way to express W is 

W^^ = TltS{k, + k2 - k3 - k4)(5(|ki|2 + |k2p - |k3p - |k4n, (3) 

where (5-functions assure conservation of the total momentum and the total kinetic energy 
(which is proportional to |kp) of incoming and outgoing particles. The collision probability, 
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expressed by = r(ki, k2|k3, > 0, is invariant under permutations {1,2} — )• {2,1}, 
{3,4} — 7- {4,3}, and {1,2} — t- {3,4}. In the present paper we will consider the case of 
three-dimensional rigid spheres with diameters a and mass m, for which F simply results in 
r'^2 = 2a'^/m [4j. For other interaction potentials, as Coulomb or Born approximation, refer 

to [iHiEn]. 

For the purposes of our work, we consider a homogeneous and isotropic (in physical space 
M^) system with the one-particle probability density function independent of x and its mo- 
mentum dependency coming only via the modulus k = |k|, so n(x, k, t) — )■ n{k, t). It is useful 
to express the distributions in the energy space Ui = |kjp where we use again the notation u 
for the energy in analogy with wave turbulence. Then, the particle density in w-space satisfies 
the relation J N{uj,t)dui = J n{k,t)dk or, in the other words, N{uj,t) = n{ui,t)^luj^ 
where is the solid angle. After these considerations Boltzmann equation ([T]) simplifies to 
the homogeneous isotropic Boltzmann equation (HIBE): 

ON f°° 

= J 'S'i2("'3'T-4 - nin2)5{ui + U2-UJ3- U}4:)du}2duj3duj4, (4) 
where we denote for brevity Ni = N{uJi,t) and = n{ui,t), and the functional 



= {UJ1U2UJ3UJ4) 



dki 



dui 



dko 



duo 



dk-i 



du% 



dki 



dud 



F?^5(ki + ks - ks - ki)dnidn2dn3dn^ (5) 

n 



takes into account the change of coordinates and the average over solid angles. Here- 
after, we always consider a three-dimensional gas of hard-sphere particles in a non- 
dimensional form with m = 1 and cr^ = 8 . Then the functional simply results in 
S'fl = 27rmin [y/uJi, y/uJ2, a/ws, y/uJl] (see Appendix |a] for details of the angular integration). 
The HIBE has two conserved quantities, the mass and energy densities, 

/ + 00 f + OO 

n{u,t)dk = 2tt / n{u,t)y/ujduj, 
-oo ^0 /^\ 

/+00 r+oo ^ V / 

n{u,t)k'^dk = 2n / n{u,t)u^du. 
■oo Jo 

Note that pm and pe are always constant in time for any distribution n and interaction 

potential, due to the fact that collisions are 2 — 2 and elastic. This is evident by evaluating 

their time derivatives using equation (|4]): the symmetries with respect to the integration 

indices immediately show that these quantities are zero. 
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A. Steady solutions 



1. Equilibrium in a closed system 

The HIBE Q is an integro-differential equation with no general analytic solution. It 
is easy, however, to look for steady (time independent) solutions. In closed system, i.e. 
without forcing and/or dissipation mechanisms, the only steady solution corresponds to the 
thermodynamic equilibrium described by the Maxwell-Boltzmann (MB) distribution, 

nMB{i^) = e T = Ae (7) 

where A = e~T and constants jj and T have the meaning of the chemical potential and 
the temperature respectively (we consider the natural unit system, where the Boltzmann 
constant is one). Validation is trivial by plugging ([T]) into (|4]): for any value of T, /i and the 
interaction potential 3^2, the (5-function assures that the integrand is zero. Moreover, the 
total mass density of the system is pm = A (vrT)^, the total energy density is pe = '^Att^T^, 
and any other moment of uj, due to the bi-parametric nature of the MB distribution, is a 
function of pm and pe- The H theorem states that in a closed system any out of equilibrium 
distribution with defined mass and energy densities will always relax to the MB distribution 
having same pM and pe- 

In Fig. [T] we show a numerical simulation of the HIBE with initial condition given by 
a Gaussian function centered around a particular value of energy; as it is clear from the 
figure, the initial condition relaxes to the MB distribution. The numerical algorithm used 



to perform this simple example will be discussed in Section |IVB[ We can observe that the 
initial condition evolves reaching an equilibrium MB distribution: the exponential behavior 
become evident by observing the inset where we plot it lin-log plot scale. Moreover by fitting 
the results with the MB function we can find the thermodynamic quantities A and T: those 
correspond exactly to ones expected knowing initial mass and energy densities (note that 
now integrals ^ are evaluated from to a finite value of lu due to numerical finiteness of 
w-space). 
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FIG. 1: Numerical computation of HIBE with an initial Gaussian shaped distribution (continuos 
black line) : intermediate states are shown with gray lines and final steady distribution with dashed 
black line. The latter has a MB behavior ([T]) with fitted parameters A and T printed in figure. 
The inset shows the same plot in lin-log scale. 

2. Nonequilibrium steady states 

Now, what can we expect in an open system driven by external forcing and dissipation 
mechanisms? We will answer this question keeping in mind the main results of the wave 
turbulence theory. Part of this theory is dedicated to study steady solutions to kinetic 
equations in the power-law form, n{uj) ~ u~^, where the constant x assumes different values 
depending on the considered wave system. It is sometimes possible to find the so-called 
Kolmogorov-Zakharov (KZ) solutions nxzi'^) ~ which correspond to constant fluxes 
of conserved quantities through scales. The KZ distribution always appears in a range of 
scales, known as inertial range, between the forcing and dissipation were the source and sink 
are located. 

As already mentioned, the HIBE conserves the number of particles and the energy, and 
so one could expect to observe two turbulent KZ cascades. The KZ exponent x can be 
evaluated by applying the standard Zakharov transformations [5j, by dimensional analysis 
[21] . or by using the method (equivalent to Zakharov transformation) proposed by Balk 
[22] • We have chosen the last one and the complete analytical calculations are presented in 
Appendix [Bj The KZ exponents depend on the scaling behavior of the scattering term 
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and on the dimension d of tlie particle system. For the particular case of three-dimensional 
hard spheres we have 

7 

constant particle flux 77 =^ riKzi^) ~ w"*, 
constant energy flux e =^ riKzij^) ^ ^ 

The simplest way to mimic an open system where steady nonequilibrium distributions of 
the form of turbulent KZ solutions can be establish is to consider a forced-dissipated HIBE 

dN 

— = hM{u:i,t) + F{u^) - D{uji)N{ui). (9) 

The forcing F is constant in time and very narrow near a particular energy value Uf-. with this 
choice the incoming fluxes of particles rj and energy e roughly satisfy relation e = Ufi]. The 
dissipation term D is implemented as a filter which removes, at each iteration time, energy 
and particles outside of the domain u G (wmm, i^max)- Further details on the numerical 



scheme are explained in Section IV B What happens if we try to solve numerically such 
forced/damped integro-differential equation? 

In Fig. |2] and Fig. [3] we plot the nonequilibrium steady states obtained with numerical 
simulations of the HIBE with forcing and dissipation; the initial conditions are characterized 
by n{uj, t = 0) = . The parameters in the simulations are Wmin = 5, Umax = 195, the forcing 
rate F = 10~^. In Fig. [2] forcing is located at loj = 22 and in Fig. |3]at = 182. No 
power-law distributions, and so no KZ solutions ([s]), are observed (note that both plots are 
in lin-log scales), but instead one can see weakly perturbed exponential curves. We can 
attempt to measure the quantities T and A in ([T]) by fitting our numerical curves; however, 
those are not perfect straight lines (in the lin-log plot) and left and right branches with 
respect to forcing scale may give different results. For such reason we will denote by 
the quantities evaluate on the left brach and with (■)/? the right ones. 

Another example we analyze is the case where we fix the forcing and dissipative scales 
and change the forcing rate. Numerical results for final steady states evaluated for three 
different forcing amplitudes, F = 10~^, 10~^, 10~^, are presented in Fig. |4j The effect of 
increasing the amplitude F results in an upward shift of the curves. Therefore, qualitatively, 
the temperature appears to be the same for each value of the flux. The only difference is 
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Al=1 .6e-04, tL=4.3e+01 
AR=2.0e-04, Tp|=2.8e+01 




FIG. 2: An example of HIBE ([9]) steady state shown with black Ime in lin-log scale. Simulation 
parameters are: Wmin = 5, (Jj = 22, Wmax = 195 and F = 10~^, and uJcutoff = 200. The dashed 
and point/dashed lines are left and right branch best fits obtained with the MB distribution ([7|): 
fitting parameters are reported in label. 




50 100 150 200 

(0 

FIG. 3: An example of HIBE ^ steady state shown with black line in lin-log scale. Simulation 
parameters are: Wmin = 5, = 182, t^max = 195 and F = 10~^, and oJcutoff = 200. The dashed 
and point /dashed lines are left and right branch best fits obtained with the MB distribution ([T]): 
fitting parameters are reported in label. 
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FIG. 4: Steady states of HIBE ([9]) in lin-log scale obtained for different values of the forcing rate 
F. Parameters are u^cutoff = 200, Wmin = 5, w/ = 21 and Wmax = 95. 




10^ 10^ 10^ 
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FIG. 5: Total energy densities PE^t) in function of time for different forcing rates F. For the system 
parameters refer to ones in Fig. |4] 

the speed at which the system, initially empty, reaches its steady state. In fig. [5] we show 
the energy density evolution (same line styles corresponds to same systems). 

After these preliminary numerical results, a lot of questions can be posed. Why no 
KZ constant flux solutions are observed but just small deviations from MB distributions? 
What happens when forcing or dissipation scales are changed? What is in general the 



10 



relation between physical quantities such as fluxes, forcing and dissipation scales and the 
MB parameters? The aim of this manuscript is to provide explanations to such phenomena 
and answer these questions. 

B. Locality of interactions 

For the KZ spectra to be valid mathematical (and therefore physically relevant) solutions, 
it is necessary that they satisfy the locality condition. A spectrum is local when the collision 
integral converges. In other words, non-locality means that the collision integral is not 
weighted scale by scale but most of the contributions come from the limits of integration 
corresponding to the ends of the inertial range. Physically, the non-locality is in contradiction 
with the assumption that the flux of the relevant conserved quantity in the inertial range is 
carried only by the nearest scales. Mathematically, locality guaranties that the KZ spectrum 
is a valid solution in an infinite inertial range, which is not guarantied a priori because 
Zakharov transformation is not an identity transformation and could, therefore, lead to 
spurious solutions. 

For the HIBE case, locality depends on the particular interaction potential, which affects 
the scaling of F'^g^ the dimensionality of the system - for detailed calculations see 

AppendixjB] Locality is not always found for both KZ solutions: for example for the Coulomb 
potential only the energy cascade is local, as shown in [T^. In the case of three-dimensional 
hard spheres considered in the present work, the criterion of locality is never satisfied for 
any of the two KZ solutions, which means that these solutions are un-physical and irrelevant 
in this model. 

C. The flux directions 

Besides locality, another important requirement for establishment of the KZ spectra is 
the correctness of the flux directions for the respective conserved quantities. In a system 
where two quantities are conserved, the following Fj0rtoft-type argument is used to establish 
which quantity must have a direct or an inverse cascade. 



11 



1. The Fj0rtoft argument 



Consider an open system where forcing scale is widely separated from a low-w dissipa- 
tion scale cjmin and a high-ci; dissipation frequency cjmax, thus ^ ^ c^max- Because 
the energy density in the w-space is different from the particle density by factor u, the 
forcing rate of the energy e is related with the forcing rate of the particles r/ as e ~ ^fV- 
Suppose that some energy is dissipated at the low scale Wmin at a rate comparable with 
the forcing rate e. But then the particles would have to be dissipated at this scale at the 
rate proportional to e/wmin ~ V'^f/'^mm ^ ^, which is impossible in steady state because 
the dissipation cannot exceed the forcing. Thus we conclude that in the steady state the 
energy must dissipate only at Wmax- By a symmetric contradiction argument one can eas- 
ily show that the only place where the particles can be dissipated in such systems is Wmin- 
This means that energy must have a direct cascade (positive flux direction) and particles an 
inverse cascade (negative flux direction). 

2. Flux directions in the HIBE 

It has been proved in [18], see also Appendix |B] for details, that fluxes of the KZ solutions 
for all types of the interaction coefficient F have always the wrong directions with respect to 
the Fj0rtoft argument requirements (in the case x > 0). An alternative way for finding the 
sign of the fluxes is considering them for general (not necessarily steady) power-law spectra 
n{u,t) ~ and plotting them as functions of x for a fixed {uj,t), see Fig. [6} Three 
exponents x correspond to steady solutions of HIBE: the particle equipartition Xeq = 0, the 
KZ particle cascade and the KZ energy cascade x^. As shown in Appendix [B| we know 
that e = on the particle cascade, 77 = on the energy cascade, whereas in the equipartition 
both fluxes are zero, i.e. e = rj = 0. We also know that for large negative x (large positive 
slope) both fluxes must be negative, as such a steep unsteady spectrum would evolve to 
become less steep, toward equipartition. Note that always x^ < x^, when x > 0. Now we 
can sketch the particle and energy fluxes as function of the exponent x as it is done in Fig. 
|6} From this sketch, it can be easily understood that whenever the condition x^ < x^ is 
valid, the particle flux will be positive and the energy flux will be negative, contradicting the 
Fj0rtoft argument. This means that one cannot match these formal KZ solutions, obtained 
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FIG. 6: Energy and particle fluxes on power-law solutions n{uj,t) ~ w ^ as functions of x for the 
three-dimensional hard sphere model. 

for an infinite inertial range, to any physical forcing or dissipation at the ends of a large 
(but finite) inertial range. 

What is then happening when fluxes have wrong direction? It has been observed in 
optical wave turbulence [I2] that the pure KZ spectra are not established in these cases 
and one has to expect a mixed solution where both a flux and a thermal components are 
present. Such mixed states are quite common for turbulent systems of different kinds, 
including strong Navier-Stokes turbulence and have been named warm cascades [191. Such 
cascades were obtained within the Leith model (which belongs to the class of the differential 
approximation models) as exact analytical solutions. 



III. DIFFERENTIAL APPROXIMATION MODEL 



Numerical integration of the Boltzmann collision integral is very challenging because the 
number of degrees of freedom grows as a polynomial. A great simplification comes from 
the isotopic assumption, which reduces the degrees of freedom from to N"^ {N is the 
number of points needed to describe the distribution). However spanning a large number 
of momentum scales is still difficult. For those reasons, some approximations to the kinetic 
equations were proposed in order to increase the range of modeled scales, see for example 
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A great simplification is to replace the collision integral operator of the kinetic equation by 
a nonlinear differential operator which mimics the basic scalings of the original one and yields 
the same steady solutions. The HIBE then results in a nonlinear partial differential equation 
called the differential approximation model (DAM). Such models have been proposed to 
simulate turbulence in different research fields: for example in water waves [21], in nonlinear 
optics [12], in strong Navier-Stokes turbulence [251 [26], iii Kelvin quantum turbulence [27], 
in astrophysics (Kompaneets equation) [28], in semiconductors [29]. Replacing the integral 
operator by a differential one amounts to assuming locality of the scale interactions, which 
means the relevant distributions must be local for DAM to have a good predictive power. 
We mentioned in Section |TT] that for hard sphere Boltzmann equation the pure KZ spectra 
are non-local and so no DAM would be advisable. However, we observed in some examples 
(Fig. [2] and Fig. |3]) that the relevant solutions in this case are not pure KZ spectra but 
distributions which are close to MB, warm cascades, which appear to be local. Thus, we use 
the DAM for describing this system, after which we will validate our results by computing 
the full HIBE. 

For the dual cascade systems, such as gravity water waves [30j, nonlinear Schrodinger 
equation ^12j, two-dimensional hydrodynamic turbulence [26], Kelvin waves [27] or HIBE 
considered here, DAM has always the form of a dual conservation law, 

dtN{u,t) = d^^R[nicu,t)], (10) 

where i? is a nonlinear second-order differential term whose details depend on the particular 
model. This equation can be written as a continuity equation for the particle invariant, 

dtN{u,t) + d^r]{u,t) = 0, 

with the particle flux 

7]iuj,t) = -d^R[n{u,t)]. (11) 



Moreover, equation (10) can be written as a continuity equation for the energy 

dt[N{co,t)co]+d^e{LO,t) = 0, 

with the energy flux 

e{u, t) = R [n{u, t)] - ud^R [n{u, t)] . (12) 
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We are now able to find the functional R by requiring it to yield the MB distribution ([T]) 



and the KZ spectra l\8n as steady state solutions of DAM (10). These constraints lead to 



R[n{uj,t)] = -Sco^ n^{uj,t)d^^ log n{uj,t), (13) 

where 5 is a constant. A formal derivation starting from the kinetic equation can be obtained 
following [121 122] • It is trivial to verify by substitution that KZ solutions ([s]) correspond 
to constant fluxes through scales. Namely, the KZ particle cascade has a constant particle 
flux and zero energy flux while the KZ energy cascade viceversa. Let us again consider the 
the flux directions on the KZ distributions, but now using DAM. Substituting power-law 



spectra n = cu ^ into (13), equations (11) and (12) yield 



7] = c^Sx(9/2-2x) Jl''-''^ 

(14) 

e = c2^x(ll/2-2x)w9/2-2x_ 

By plotting r/ and e as functions of the exponent x at fixed w, we arrive again at Fig. |6| 
Note that it is by using DAM such plot was obtained. Once again we note that the particle 
and the energy fluxes on the respective KZ solutions [x = 7/4 and x = 9/4) have wrong 
directions with respect to the Fj0rtoft argument. 

The beauty of the DAMs is the possibility to solve numerically the system for wide 
frequency ranges and, therefore, to find clear scalings. In particular, such models are very 
efficient for finding constant steady flux solutions because they become simple ordinary 
differential equations (ODEs). In the following we will present some analytical results for 
such steady states. 



A. Constant energy flux: direct cascade 



We will now find an ODE that describes a constant direct energy cascade e with no flux 
of particles, which we call ODE-e. According to Fj0rtoft argument, this implies a large 



direct-cascade inertial range. Putting 77 = in (11) and (12), we have 



constant energy flux 



e = R{ll!, t) = const. 



(15) 
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Using (13), we arrive at the following Cauchy problem 



n{ujo) = no, (16) 

where we have chosen the boundary conditions fixing the values of the distribution and its 
derivative at the same point uq (e.g. at the forcing scale) for ease of numerical solution. 

If we solve numerically in w-forward the ODE-e for different values of the energy flux we 
flnd curves presented in Fig. [7j Here we do not want to discuss the details (it will be done 



widely in Section IV), but just remark that the solutions follow the MB distribution and 
suddenly change behavior going very fast to a zero value of the distribution. We will call 
this rapid change a front solution. 

1. Compact front behavior 



It is possible to flnd a front solution for the equation ( 15 ) describing the behavior near 



the dissipation scale. Let us seek for a front solution which in the vicinity of a certain point 



i^max behaves like n{uj) = B (wmax — ■ If we plug this expression into (15) and take the 
limit u — )■ Wmax "we flnd that to satisfy this equation in the leading order in (cjmax — ^) we 
must have 

a = 1 I I 

= A/ ^ IS/2 (^max-C^)- (17) 



R — / t ^ ' \ Q 13/2 

-D — \ Q 13/2 V *J t'-'max 

y -J t'Jmax 

Thus, the front solution is linear in the vicinity of Wmax with a slope depending on the 
dissipation scale cUmax and the value of the energy flux e. Note that the compact front 
behavior at the dissipation scale is typical for DAM. We will soon discover that Wmax is a 
very useful physical parameter which allows us to flnd a link between the temperature, the 
chemical potential and the energy flux in the forced-dissipated system. 

2. Kats-Kontorovich correction 

Lets summarize our preliminary observations. We expect a warm cascade, that is a 
distribution which contains both the flux and the thermal components. We have also found 
that the solution has a compact front which arrests the cascade at the dissipation scale 
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Wmax- We will now assume (verifying it later) that in the most of the inertial range the 
warm cascade solution is close to the thermodynamic MB distribution and the correction 
due to finite flux is small. We then perform a qualitative matching of the fiux-corrected 
MB distribution to the compact front, and thereby obtain a relation between Wmax, T and 
A in 0. To find the warm cascade solution in the inertial range, we consider the Kats- 
Kontorovich (KK) correction to the Maxwell-Boltzmann distribution: 

n{u) =nMB{l + n) = {l + n)Ae~^, (18) 



where n is small, n ^ 1. By plugging this solution into (15) and linearizing in n we end up 
with the following ODE-e for the correction 

euj'^A'^eT^ = -SdujujU. (19) 



3. Matching 

We will now match the KK correction to the front solution. The basic idea is to force 
the KK solution to satisfy the n(co'inax) = and to have at Wmax the same slope as the front 
solution. Detailed calculation is presented in Appendix [C] The prediction results in: 

e = ScoLA'e-'^. (20) 

This relation is very important because it gives an analytical relation between the thermo- 
dynamic quantities T and A in terms of the energy flux e and the dissipation scale Wmax- 
However we note that our matching is only qualitative, because the KK correction is sup- 



posed to be small which is not the case near the front. Thus, the relation (20) is approximate 
and we do not expect it to hold precisely. 



4- Alternative approach to find Wmax 

Another simple way to find a prediction for the value of Wmax is the following. As we 
expect to observe a warm cascade, we can ask what will be the range where the thermal 
component will dominate the dynamics. We can simply assume that in most of the inertial 
range we will have a distribution n ~ umb- Note that the MB distribution always has a 
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positive concavity, dujui^ > 0. On the other hand, we note that our ODE-e can be re-written 

as 



d,.„„n 



1 

n 



13 

b 00 2 



(21) 



from which it is clear that dujujU may change sign. The point at which dujujU = can be 
considered as s boundary separating the MB range (with negligible flux correction) and the 
front solution (with large flux correction). This boundary can be estimated by a simple 



substitution of the MB distribution to the r.h.s. of (21), which gives 



S u 2 e T 

2^2 



g,{uj,A,T). 



(22) 



As this relation contains the exponential factor which decays very fast (for Wmax ^ T, see 
Appendix [C]), it is natural to think that the range at which e becomes important appears 
very sharply and is very near to the point Wmax- Thus we arrive at the following estimate. 



A-^ S CUmax e T 
2^2 



(23) 



B. Constant particle flux: inverse cascade 

In analogy of what has been done for the direct cascade, we now look for predictions 
in the inverse particle cascade r] with no flux of energy. The ODE-r/ that describes such a 



cascade is simple to obtain: by integrating equation (11) once and putting e = in (12), we 
have: 

. 1 n R(uJ,t) 

constant particles nux =^ r] = = const. 



UJ 



(24) 



This yields the following Cauchy problem. 



Tj = Suj^n^{u)d^^ logra(w), 
n{ujo) = no. 



(25) 



This problem is most naturally solved backwards in the w-space, as we are interested in 
the inverse cascade. We seek for a solution having a particle flux going from high to low 
frequencies, i.e. rj < and for convenience we will make the substitution rj — t- — 1?7| in 



equation (25). The Cauchy problem (25) is very similar to (16) with the only difference in 



the cj-scaling. Thus we will use the same approach for studying it. 
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1. Compact front behavior 



Let us find a front solution for the equation (24). We now expect the front to be on 



the left edge of the (inverse cascade) inertial range, i.e. in the vicinity of a certain point 



i^min < ^f- By plugging n{uj) = B {u — Wmin)'^ expression into (24) and taking the limit 
UJ — )• Wmin, in the leading order in {co — Umm) we have 



a = 1 



= J ^ '^11/2 (^ -^min). (26) 



Thus, the front solution for the inverse particle cascade is also linear in the vicinity of cjmin, 
with a slope depending on cjmin and the value of the particle flux t]. 

2. Kats-Kontorovich correction 

As previously supposed for the direct energy cascade, we expect in the most of the inverse- 
cascade range a corrected thermodynamic spectrum and a front solution behavior at the left 



end of this range. Let us evaluate the Kats-Kontorovich correction (18), and after that 



match it to the front solution. By plugging the expression (18) into (24) and linearizing in 
n we obtain the following ODE-rj for the correction. 



\t]\uj-^ A-^e^ = -Sd^^h. (27) 

3. Matching 

Again, we want to match the KK correction to the front solution. The idea is very similar 
to the previously used for the direct cascade, except for the fact that now the limit taken is 
'^min ^ T] for details refer to Appendix |D} This results with the following condition on the 
flux, 

\V\ = S (I) A'uL- (28) 



19 



4- Alternative estimate of uJj 



Again, we can obtain an alternative estimate for predicting the range of the warm cascade. 
Let us rewrite the ODE-rj as 

1 



n 



D a; 2 



(29) 



Keeping in mind that the MB distribution is always characterized by a positive concavity, 
i.e. d^uj^ > 0, and considering the hypothesis d^jfi ^ BujUmb we find 

^ 9 ^ 11 2cj 
A S ( I J 2 P T 

\V\= =9MA,T). (30) 

Similarly to what we have done for the inverse cascade, we now can suggest that the change 
of concavity occurs near cumm- This results in 

I-,! = ^"7/ ' . (31) 

However, we do not expect a good prediction as before because in this case the exponential 
term is not a rapidly varying function near Wmin- 

C. Double cascade 

We have now all tools to study the double cascade process. Let us force at w/, dissipate 
at Wmax and u^m, and consider the case cjmin ^ ^ Wmax- If the forcing range is narrow. 



the simple relation e = rjUf holds for the fluxes. Using this relation, and combining (20) 



and (28), we can estimate T and A in the system: 

2w 



max 




+ In - 2 In ^ 



O 7/2 ' 
'^min 



9 ' 

(32) 

^ /. I \n\ 

and, therefore, the chemical potential 

, = T(^^l„^ + l,>-J. (33) 

Note that the temperature appears to be independent of the fluxes and is completely con- 
trolled by the forcing and the dissipation scales. This means that increasing the forcing 
strength without moving Uf simply adds more particles into the system with the energy per 
particle remaining the same. 
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IV. NUMERICAL RESULTS 



In this Section we present the numerical results obtained by using the DAM and by 
integrating, at lower resolution, the HIBE. Our aim is to compare results for the warm 
cascade solutions of DAM, which has been devised as a local approximation of the integral 
collision operator, with direct numerical simulation of the full integro-differential equation 




A. DAM resutls 



We will first present some numerical experiments on integration of the Cauchy problems 



(16) and (25) in which we take for simplicity S = I. Note that all numerical simulations 
can be performed without any loss of generality starting with a particular value ujq because 
of re-scaling properties described in Appendix |E} 



1. Constant direct energy cascade 

In Fig. [T] we show the results obtained by integrating equation (16) with loq = 3.5 for 
different constant energy fluxes e. As initial conditions, we choose the values of the spectrum 
no and its slope tiq from the MB distribution having A = 1 and T = 1. The solutions follow 
the thermodynamic solution (shown as a continuous line) until they rapidly deviate and reach 
the front in the vicinity of particular values of Wmax- This numerical experiment exhibits two 
important facts always observed in simulations performed with different initial conditions: 
the presence of a long transient in which the flux correction is negligible with respect to the 
thermodynamic MB distribution and the presence of a particular value Wmax at which n{u) 



goes to zero. A lin-log plot of the function g^ico, 1, 1), see equation (22), is shown in Fig. 
Sj intersection of this curve with horizontal lines at e = 1, e = 10^^ and e = 10~^ marks 
the predicted cut-off frequencies for the respective flux values. Agreement with the behavior 



in Fig. [7| is evident: the values of Wmax obtained with equation (22) and Fig. |8j coincide 
with the observed values in Fig. [t] within 5%. Note that the peak of ge{uj, 1, 1, ) is around 
u = 3.5: this is why we set this value as initial condition cuq- 
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10 12 14 16 18 20 
0) 



FIG. 7: DAM simulations of (16) for different constant energy flux e starting with the same initial 



condition at ujq = 3.5 given by the MB distribution UMsi^) = t where A = 1 and T = 1 



(continuos line). 




In Fig. [9]we present the results for a particular case with flux e = 1. We can appreciate the 
presence of warm cascade and the front solution near cjmax- The linear behavior of the front is 
evident in the zoom near Wmax showed in the inset. Numerically we are able to measure Wmax 



and so evaluate B from equation (17). The theoretical prediction agrees with the measured 
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FIG. 9: DAM simulation of ODE-e (16) with constant energy flux e = 1 (dashed line). The initial 



conditions in = 3.5 are set by the MB distribution with T = 1 and A = 1 (continuos line). The 
inset shows a zoom of numerical n{uj) (dots) in the vicinity of the point Wmax where a linear fit is 
shown by continuos line. 



0.997%. The error is evaluated as B„ 



\B„ 



Best I / B 

n 



slope with the error B^r 
where Bmeas is the measured linear coefficient and Best is the one taken form relation (|17|). 



In all other simulations performed with different values of e or different initial conditions, 
Berr is always within 5%. 



We now check numerically the validity of the matching prediction ( 20 ) by taking different 



initial condition UMBi^o = 3.5) varying T and keeping A = 1 and e = 1: results are plotted 



in Fig. 10 It is evident from the figure that the predicted temperature (continuous black 



line) evaluated from relation (20) is an overestimation of the numerical results (dots) and 



the error is around 10%. Finally prediction for the alternative temperature relation (23) is 



plotted with gray dashed line: it appears to give a better estimation than relation (20) 



2. Constant particle cascade 



We now investigate the inverse particle cascade by solving Cauchy problem (25) going 



w-backward. In Fig. [TT]we show numerical results obtained by taking initial conditions at uq 
from MB distribution nMsi^) = Ae~^ with T = 1, A = 1. As in the case of constant energy 
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FIG. 10: Checking of predictions in DAM constant energy flux cascade (16): the points represent 
the temperature of the initial condition T with respect to the measured ti^max- Solid line is the 



matching relation (20) while dashed one is obtained from (23) 



flux, here the warm cascade range is wider for smaller flux values. We also observe fronts in 



vicinities of cutoff points Wmin- In Fig- 12 we show the function 1, 1) which represents 



the prediction of the thermodynamic range (30). Qualitative front values of results in Fig. 



11 show poor agreement with this naive estimation. 



The front solution is analysed in detail in Fig. 13 where we choose the particular case 
with ?7 = —1. The linear behavior is demonstrated in the inset. Moreover a numerical 



estimation of aJmin lets us evaluate -B, see equation (26). The error B^^t is presented in the 
figure; for all other simulations we have performed i?err remained within 4%. 

Finally we check KK matching prediction for the thermodynamic quantity A with re- 



spect to Wmin presented in equation (28): results are showed in Fig. 14 In this case the 



KK analytical prediction (continuous line) underestimates the numerical data while the 
estimation (31) is completely out of range (dashed line). However the scaling A ~ w 
KK prediction tends to be reached for small values of Wmin, where cjmin ^ T . 



W of 
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FIG. 11: DAM simulations of (25) for different constant inverse particle flux t] starting with the 



same initial condition at ujq = 3.5 given by the Boltzmann distribution nMB{'^) = t where 
A = \ and T = 1 (plotted with continuos line). 
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FIG. 12: Plot of the function 5r^(cL!, 1, 1), see equation (30) 



3. Double cascade 



An example of double cascade is presented in Fig. 15 where we set the forcing at = 
ojq = 3.5. We show here three cases where the particle fluxes are respectively 1] = — 1, 
1] = — 10~^ and r] = —10""^. Measuring cUmin and Wmax for each case we are able to estimate the 
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FIG. 13: DAM simulation of (25) with rj = —1 (dashed Hne) starting with initial condition at 
ujQ = 3.5 given by the MB distribution nMsi^) = Ae~T where A = 1 and T = 1 (plotted with 
continuos line). Inset: lin-lin scale zoom in the vicinity t^min (dots) where the best linear fit is 
presented with a continuous line. 



temperature Test from prediction (|32j). Results do not agree with the expected temperature 
(the initial conditions set it at T = 1) but they approach this value for bigger ranges, 
i.e. when the condition cjmin ^ cu/ ^ i^max is better satisfied (see for example the case 
ri = -10"^). 



B. HIBE results 



We now to present results of the direct simulation of HIBE with the full Boltzmann colli- 
sion integral and compare them with predictions obtained by DAM. As we have mentioned 
above, the evaluation of ^ is numerically challenging and it is nowadays practically impos- 
sible to simulate such wide cu-space ranges as we have done using the DAM. In the present 
work, we will always use a low resolution of 101 points by considering u G [0,Ucutoff] and 
taking a uniform distribution with Au = -^j^^- We have checked that the numerical solu- 
tions are mesh independent by taking a finer mesh, 201 points, and comparing the solution 
of one critical case. 

The 5-function in ^ defines a resonant manifold over which the integrand need to be 
evaluated; numerically it is a set of discrete resonant conditions JH = {ui, U2,oj3,oj4} / uji + 
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FIG. 14: Checking predictions in DAM constant inverse flux cascade (25): the points represent the 



thermodynamic quantity A with respect to the measured Wmin- Continuos hne is the KK matching 



prediction given in equation (28) while dashed one is obtained from (31) 



U2 = UJ3 + which can be pre-computed. Note that the dissipation at high wave numbers is 
chosen to satisfy Wmax < ^cutoffl'^ in order to prevent ultraviolet bottleneck effects. The time 
evolution is performed by using the Euler scheme. Further details on numerical methods for 
solving the HIBE and a simple code can be found in |31j . 



1. Direct cascade study 

We first analyze the direct energy cascade by putting the forcing scale near the low- 
u dissipation scale in order to have a wider direct inertial range. Numerical results for 
these final steady states were previously presented as examples in Fig. |2] and Fig. |4} 
We concentrate now only on the last one: here we kept fixed Wmin = 5, = 21 and 
i^max = 95 and varied the forcing coefficient, i.e. the fluxes r] and e. We were claiming 
that the temperature of the systems is the same because qualitatively the distributions 
have identical slopes. Moreover we observed in all the examples that left and right branch 
chemical potentials and temperatures can be defined by the forcing scale. 

With these previous DAM results in mind we have measured A and T in three examples 
presented in Fig. |4j the results are shown in Fig. 16 and are compared to analytical 
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FIG. 15: DAM double cascade simulations of equations (16) and (25) for three different values of 



the particle flux rj (and consequently of the energy flux e = riujf). The initial condition are taken 
at Ljjf = 3.5 from the MB distribution with T = 1 and A = 1. Measuring oomin and uJraa.^ in each 



case we estimate of the temperature Tg^t from prediction ( 32 ) 



predictions (32). As expected the quantity A ~ ^ but the line (in log-log plot) is shifted 
with respect to the interval between A^ and A/j, represented respectively with filled and 
empty circles. However, the theoretical prediction is much closer to A^, which is natural 
because the right inertial interval is wider than the left one. In fact, the agreement of 
Ar with the theory is quite good considering the presence of the undefined constant 5* in 



the theoretical prediction. The temperature is shown in Fig. 17 even though Tl and Tr 
are different they both appear to be forcing independent, as predicted. The temperature 



evaluated from relation (32): temperature (dashed line) stands in between of these values, 
and closer to Tr, which, again, is natural because the right inertial interval is wider. 

We have also analyzed sensitivity of the temperature to varying the high-w dissipation 



range and results are presented in Fig. 18 Keeping the forcing constant and changing the 
value of Umax the system reaches steady states characterized by different temperatures Tl 



(filled circles) and Tr (empty circles). The prediction (32), shown by the continuous line, is 
in between of the two temperatures and is closer to Tr - again due to the wider right range. 
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FIG. 16: Results for the fitted values of the thermodynamic quantity A = e t plotted against the 
forcing levels, as obtained in the simulations shown in Fig. [4| the values of are shown by filled 



circles while - by empty ones, the blue line refers to prediction (32) with 5 = 1 
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FIG. 17: Measured temperatures Tl (filled triangles) and Tr (empty ones) obtained in the simu- 
lations shown in Fig. [4| The dashed line is the analytical prediction ( |32| ). 

2. Inverse cascade study 



Finally, we have performed some simulations putting the forcing scale near the dissipation 
at high w's in order to study the inverse cascade process. In this case too, as reported in Fig. 
[3} we observe two different values of thermodynamic quantities on the left and on the right 
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FIG. 18: Temperature in different steady state keeping the forcing constant and varying the dissi- 
pation scale ujma.x- big empty circles correspond to the temperature Tr on the right of the forcing 



scale, whereas small filled circles to the left side, T^. The continuous line is the prediction (32) 
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FIG. 19: Temperature in different steady state keeping the forcing constant and varying the dissi- 
pation scale oj^in- big empty circles correspond to the temperature Tr on the right of the forcing 



scale, whereas small filled circles to the left side, Tl. The continuous line is the prediction (32). 



from the forcing. Here we are able to study the scaling of the thermodynamic quantities 
T and A with respect to changes of the small-w dissipation scale Wmm- Results for T are 



shown in Fig. 19 and for A in Fig. 20, with the "left" quantities shown by filled circles 
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FIG. 20: Thermodynamic amplitude A for different ujmm and fixed forcing in log-log scales in 
numerical simulations of HIBE. Filled circles corresponds to while empty circles - Ar. The 



continuos line is formula ( 32 ) with S = 1. 



and the "right" ones by empty circles. There is a reasonably good agreement of T with the 



prediction (32) for small Wmin- This is natural because smaller Wmin corresponds to larger 
inverse cascade inertial range and also because the prediction is valid when Wmin ^ T. 



On the other hand, for A the prediction (32) is in better agreement with the data at large 
cUmin with cjmin ^ oj f . This is duc to two possible reasons. First, we underline that the 
agreement can be made more suitable since the analytical prediction contains the undefined 
order-one parameter S which could be adjusted to better fit the numerical results. Second, 



the particle flux which defines A in relation (32) can be smaller due to finite range effects. 



Indeed, following [32], the ratio of the leftward particle flux to the total particle production 
rate is estimated in 

TjL = V ^ — —■ (34) 

This equation, in addition to other corresponding to rightward fluxes in the cited paper, 
states that for the particle flux to be mostly to the left inertial ranges in both directions 
must be large (note that this is also the condition of validity of the Fj0rtoft argument). 



Fig. 21 shows the behavior of the normalised left measured flux rjL/r] (empty triangles) with 
respect to Umin- We can clearly see that the measured particle flux is indeed much smaller 
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FIG. 21: The plot shows in log-hnear scale the ratio of the measured inverse particle flux over the 
total one rji/rj (empty triangles) with respect to Wmin; the continuous line correspond to the total 



flux while the dashed one follows the finite range correction ( 34 ) . 



than the one imposed by the forcing term (continuos line), around one third of it. This is in 
quite good agreement with the finite range prediction ( [34) ) plotted with dashed line. Similar 
reasoning can be made for corrections on Adr/I) in the case of direct cascade example in Fig. 

m 



V. CONCLUSIONS 



In the present paper we investigated stationary turbulent states in the isotropic Boltz- 
mann kinetic equation for hard spheres. This was done by looking for steady nonequilibrium 
states in open systems, that is when forcing and dissipation mechanisms are present. Analo- 
gies with similar results of wave turbulence theory suggest the manifestation of a warm 
cascade, i.e. a constant direct flux of energy and inverse flux of particles on background 
of thermodynamic Maxwell-Boltzmann distribution. This is a consequence of wrong flux 
directions in KZ solutions with respect to the Fj0rtoft argument. 

We have built an ad-hoc differential approximation model to easily simulate the cascade 
processes. Indeed, this simplification allowed us to reach a wide range of scales inaccessible by 
solving the isotropic Boltzmann kinetic equation directly. Simulations show the presence of 
a warm cascade with approximately the MB shape followed by a sharp front for both energy 
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and particle cascades. We have physically interpreted o^min and o^max as intrinsic dissipation 
scales at low and high w's which are necessary to estabhsh the steady state. Moreover, 
we have found analytical predictions relating the particle and energy fluxes, forcing and 
dissipations scales to the thermodynamic quantities of the system. In particular we have 
shown that the temperature is independent of the amplitude of the fluxes but only depends 
on the forcing and dissipation scales. 

We have then compared the theoretical predictions and the numerical results obtained 
with the differential approximation model with simulations of the complete isotropic Boltz- 
mann kinetic equation. Even though the resolution for the latter was limited by the available 
computational power, the results are comparable and in good agreement with the analytical 
predictions. In particular we have verified that the steady state is characterized by a warm 
cascade where a fitted thermodynamic Maxwell-Boltzmann distribution has been used to 
measure temperature and chemical potential of the system. We observe, in agreement with 
our analytical predictions, that the temperature is completely defined by the forcing and 
dissipation scales and does not depend on the fluxes. 

We hope that this work may open some perspectives towards understanding nonequi- 
librium steady states and their net currents (fluxes) by cross-fertilization with the weak 
turbulence theory. 
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Appendix A: Three-dimensional (5-function angular average 

The angular average of the four- wave hnear momentum conservation 5(kf2) = S(k.i + 
k2 — ka — k4) is evaluated by splitting it into two 5-functions of three particle collision. This 
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results in 



/ 5i]4l)dn,234 = 5(ki + k2 - k)5(k3 + k4 - k)dkdnu34 

= / 5(ki + k2 - k)c/fii2 / (5(k3 + k4 - k)dn 



ffc: 



1 1 .9.. 27r 



where geometrically fcmin = l^i — ^2! = 1^3 — ^4] and /cmax = ^1 + ^2 = ^3 + ^4- For details 
about the integration of three particle (5-function see Appendices in [5]. 

Appendix B: Kolmogorov-Zakharov solutions for general HIBE 

The Boltzmann collision integral Icoii is defined as 

/oo 
Til [n(x, ka, t)n(x, k4, t) - n(x, ki, t)?2(x, k2, t)] 
'OO 

X5(ki + k2 - k3 - ki)5(|ki|2 + |k2p - Iksl' - |k4nrfk234, (Bl) 

where the two (5-functions assure the conservation of the linear momentum and kinetic energy. 
In the isotropic case it is convenient to move in the energy domain = |kjp G [0, +00) and 
so the HIBE results in 



00 

/(Wi) = / S'f2(^3^4 - ^in2)5(a;i2)'^^234, (B2) 







d-1 



where /(wi) = fix Jcoii(x, wi, ^ ^ and we use for brevity = n{uji) = n(x, |kjp,t), 
and ^(co'il) = S{ui + U2 — C03 — UJ4). The functional S is 

Sfl = hu,U2Usco,y^-'{Tlt6{k, + k2 - k3 - k4))n (B3) 
lo 

and the operator states for the integration over solid angles. It is important for the 
following to estimate the homogeneity degree of S. Supposing that the coUisional kernel 
scales as r^|'^2) = ^"^^^li^ have 

Sxlu) = X^'-''^-^'^"- = x'^'-'^-'Sfl (B4) 

Moreover its behavior at the boundaries of integration is 

lim Sft-^t'^^' 
^'"""^ . (B5) 

lim ~ 

34 



if we assume that 



and 



lim (r?^5(ki + k2 - ks - k4))n ~ (B6) 

UJi—^+OO 



lim {Tlt5{k, + k2 - kg - k,))n ~ coj' (B7) 

(note that for cjj — )■ oo also another Uj must go to infinity due to the ^-function). 

In the following we will suppose that the particle distribution function follows the power- 
law distribution n{uj) = Auj~'^ and so 

/(Wi) = I 5'^3+4_i) [oJ^^UJl" - UJi" (W3 + W4 - Wl)""] 0(^3 + W4 - ^i) du^^i, (B8) 

^0 

where is the Heaviside step function. 



1. Kolmogorov-Zakharov solutions 

We will present the Kolmogorov-Zakharov solutions of the collision integral using the 
method presented by Balk in [22] . The collision integral, without any loss of generality, can 
be rewritten as 

J(u;i) = A^uj^^-^ r Sit - {^i^2U^z^a) c^f (B9) 

Jo ^2 ^2, 

where the exponent 

/i = 2z/ + l-2/3-y (BIO) 

is chosen in order to have zero as homogeneity coefficient of the integrand (excluding the 
differentials ^). If the integral converges. Balk proved that is possible to interchange the 
three integration index in the integrand with the fourth one, Ui. Thanks to the symmetric 
properties of the collision kernel we can write 

4 Jo 

x(wi +u;^-W3 -Od(a;i2) > (BH) 

UJ2 OJ3 OJa 

which clearly vanishes for /i = or /i = 1. This corresponds to the condition on the exponent 

3d -2 ^ 
^0 = i^U=o = — \- P 

yi = i'\^=i = — + p. 
Note that first KZ solution for HIBE were presented in [17j . 
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2. Convergence of the integral (locality condition) 



The locality of interactions is guaranteed by the convergence of the collision integral. We 
then investigate the possible values of v which assure the convergence around the integrand 
singularities. 

a. Limit — )• oo 

In the limit of c^a — )■ oo we can approximate (ws + — = u^'^ — utu^'^^^ (104 — ui) + 

0{u^'^~'^) at the second order. The argument in the square brackets of (B8) results in 



- W3' 



[(JJ4 " — + "lo^ ^{ui — LOi)] . (B13) 



As a consequence, when z/ > 0, the integrand for large U3 goes like and so the 

convergence condition is 

jy>d-l + Ti. (B14) 



b. Limit ws — )• 0"^ 

In the limit of ^3 — )■ 0+ we can approximate (c^s + co^ — Ui)^'^ = {u^ — Wi)^'^ — uusi^u^ 



cui) " + 0{ujI) at the second order. The argument in the square brackets of (B8) results 
in 

[...] = Lo^" [tu^" - tu^'co^iui - ui)-" + z/u;r"w3^'(w4 - tuiy-^] . (B15) 

So, when i/ > 0, the integrand for small Go's goes like and so the convergence 

condition is 

iy<^ + r2. (B16) 
Analogue condition holds for the singularity (ws + — Ui)''^ ^ 0^ . 

3. Constant fluxes 

The solutions n{u) = Au'"'-^ and n{u) = Au^"^ correspond, respectively, to constant 
flux of particle and energy. To demonstrate this fact we perform the substitution Ui = Ui^t 
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Vi 7^ 1 in the equation (Bll) which results, recalhng the homogeneity of the integrand 
function, in 



J A 



X (1 + ^2^ - ^3^ - ^4^) 5{l + 6 - 6 - ^4) C?634 = ^'"^Z ' Uifl) (B17) 



The integral U (/i) is now performed over the triangle A in the .^3 x ,^4 space satisfying the 
conditions < < 1 and ^4 > 1 — ^3, without any dependence on ui. 



a. Flux of particles 

The flux of particles is defined as 

Q{u}) = - /(wi) dui = / ^ ''dui = P . (B18) 

Jo 4 Jq 4/i 

If /i = 1 the flux is zero while in the case /i = it is indeterminate. By applying the De 
I'Hopital rule in the latter case we find 

Qico) |,=o = tL " ^^'^'^'^ (Ik) + " " ^^^^^ 

The integrand, and so the sign of the particle flux, is always negative for i/q > 0. This is 
clear by looking at the sign of every factors in the integral: all are trivially positive except 
- i^""") and In (^) . Recalling that (1 - - ^4) > and 6 = ^3 + ^4 - 1 we 

have 

< (1 - - ^4) = 6^4 - 6 - ^4 + 1 = - 6 ^ 6^4 > 6, (B20) 

which leads to In ^ and {^^^'^^^^^^ — ^2"°) — positive z/q). As a consequence 

Q{uj) > 0, that is the particle flux goes from low to high frequencies. 



b. Flux of energy 
The flux of energy is 

P{uo) = - / /(wi) uji dui = / ^dui = \ (B21) 







4 7o ' 4(1 -/i) 



37 



and is null when fi = while indeterminate in the case /i = 1. Again applying the De 
I'Hopital rule we have 

Pi^) l,=i = ^ ^ Sfll^ (^3-'^^^'^^ - er^) (66^4) (B22) 

X [6 ln(6) - ^3 H^S) - U HU)] 5(1 + 6 - 6 - ^4) ^^634 

As previously discussed, the term (^3''^^^^'^^ — ^2^*^^ < for every z/i > 0. Differently, the 
factor [^2 111(^2) — ^3 111(^3) — ^4 ln(^4)] is always positive but here the demonstration is not 
so trivial as in the previous case and for a complete discussion see [18]. So P{u}) < 0, which 
means that the energy flux goes from high to low frequencies. 



Appendix C: Matching Kats-Kontorovich to ujmax front solution 

We will here find the match between the KK correction and the front solution for the 
ODE-e. We make the hypothesis that the front occurs for Wmax ^ T and so it is reasonable 



to think that the term 00 2 in equation ( 19 ) it is slowly varying with respect to e ^ . So by 



integrating twice in a; (19) and match to the front we get the Cauchy problem 



_ 13 n^2 . o 2a; ^ ^ y % 

eu 2 —A T = -t,n + ci[uj - Wmax) + C2 



^(Wmax) = -1 (CI) 
<9u;^(c^max) = B 

where the first condition assures that niujm»^) = and the second that the front behavior is 



linear with slope B = —eiUmax found in equation (17). The integration constants are then 

13 



*;max 



O I ~"2" T2 ,1— 2 

C2 = -S +€ Wmax A ^ 6 T 

_ 13 
) 2 

2 



4 

Ci = S 2 62 Wmal A '^e T - e Wmax ^ A ^ 6 



13 13 _ _ , (C2) 



max 



T 



We will now match this solution in the regime where T ^ Wmax and the flux is negligible 
with respect to the thermodynamic solution. In this regime, by assuming the scaling relation 
e ~ UmaxA'^e , the coefficients results in 

" (C3) 

„_1 i — ^ "max 

and so the smallness of the correction reads as 

~ , , * O 2aj . \ ^ 1 1 T- t 1 "max ™ / ^ , \ 

n[uj) = = -euj"^ —A'^e^ + (w - cJmax)5"2e2a;maxA" - S. (C4) 
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Finally, considering that uj <^ Wmax, we recover and validate the relation 

e = ScoLA'e-'^. (C5) 

Appendix D: Matching Kats-Kontorovich to comm front solution 



The KK correction for that case is given by equation (27). We will consider the limit 



u <^T and so e r ~ 1. By integrating twice in u we get the Cauchy problem 

< n(a;min) = -l (d) 

1 - — 

with the initial conditions chosen in order to match with the front solution. The integration 
constants result in 



7 
' 2 

'min 63 

_ 11 

ci = -St \ri\iuj^-^ + \ri\ujJjA 



C2 = -S+ \r]\ujJ^^^A 2 



mm ' I / 1 mm 9 ^ 

7 

We assume and guess that the particles flux scales as \ri\ ~ u^^^A"^. Now, in the regime 
u ^ Wmin were the correction is negligible we have 

7 4 

n{uj) = = -|?7|u;^2 + ci{uj - Umm) + C2 ~a;ci. (D3) 

00 

So, finally, we impose that ci = to get the condition on the flux 

\v\ = s (I) A'uL- m 



Appendix E: Scaling properties of DAM 



The constant energy flux DAM (15) and the constant particles flux DAM (24) can be 
generally written as 

c = -SujPn'^{u)d^^ log n{uj) (El) 

where the exponent is respectively p = 13/2 and p = 11/2 and c is a constant that represent 
the considered flux. Lets now analyze the rescaling properties of that equation by the 
following change of variables 

c= A"c 

u = X^co (E2) 
n = X'n. 
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After some easy algebra we find that the system is invariant if 

a = {p-2)/3 + 27. (E3) 

As a consequence we can estabhsh how the thermodynamic quantities defined by the 
Maxwell-Boltzmann distribution nMsi^) = Ae~^ = e~^^ vary: the temperature T scales 
as u and so T = X^T, while the chemical potential /i scales as = X^jft. 
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